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Abstract 

This paper is about a curious phenomenon. Suppose we have a data matrix, which is the 
superposition of a low-rank component and a sparse component. Can we recover each component 
individually? We prove that under some suitable assumptions, it is possible to recover both the 
low-rank and the sparse components exactly by solving a very convenient convex program called 
Principal Component Pursuit; among all feasible decompositions, simply minimize a weighted 
combination of the nuclear norm and of the l\ norm. This suggests the possibility of a principled 
approach to robust principal component analysis since our methodology and results assert that 
one can recover the principal components of a data matrix even though a positive fraction of its 
entries are arbitrarily corrupted. This extends to the situation where a fraction of the entries 
are missing as well. We discuss an algorithm for solving this optimization problem, and present 
applications in the area of video surveillance, where our methodology allows for the detection of 
objects in a cluttered background, and in the area of face recognition, where it offers a principled 
way of removing shadows and specularities in images of faces. 

Keywords. Principal components, robustness vis-a-vis outliers, nuclear-norm minimization, 
£i-norm minimization, duality, low-rank matrices, sparsity, video surveillance. 



1 Introduction 

1.1 Motivation 

Suppose we are given a large data matrix M, and know that it may be decomposed as 

M = L + So, 

where Lq has low-rank and So is sparse; here, both components are of arbitrary magnitude. We 
do not know the low-dimensional column and row space of Lq, not even their dimension. Similarly, 
we do not know the locations of the nonzero entries of So, not even how many there are. Can we 
hope to recover the low-rank and sparse components both accurately (perhaps even exactly) and 
efficiently? 

A provably correct and scalable solution to the above problem would presumably have an impact 
on today's data-intensive scientific discovery^] The recent explosion of massive amounts of high- 

1 Data- intensive computing is advocated by Jim Gray as the fourth paradigm for scientific discovery [24]. 
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dimensional data in science, engineering, and society presents a challenge as well as an opportunity 
to many areas such as image, video, multimedia processing, web relevancy data analysis, search, 
biomedical imaging and bioinformatics. In such application domains, data now routinely lie in 
thousands or even billions of dimensions, with a number of samples sometimes of the same order 
of magnitude. 

To alleviate the curse of dimensionality and scalej^jwe must leverage on the fact that such data 
have low intrinsic dimensionality, e.g. that they lie on some low-dimensional subspace [15], are 
sparse in some basis [13], or lie on some low-dimensional manifold [4,46]. Perhaps the simplest and 
most useful assumption is that the data all lie near some low-dimensional subspace. More precisely, 
this says that if we stack all the data points as column vectors of a matrix M, the matrix should 
have (approximately) low-rank: mathematically, 

M = L + N , 

where Lq has low-rank and Nq is a small perturbation matrix. Classical Principal Component 
Analysis (PCA) [15, 25, 27] seeks the best (in an £ 2 sense) rank-fc estimate of Lq by solving 

minimize \\M — L\\ 
subject to rank(L) < k. 

(Throughout the paper, ||M|| denotes the 2-norm; that is, the largest singular value of M.) This 
problem can be efficiently solved via the singular value decomposition (SVD) and enjoys a number 
of optimality properties when the noise Nq is small and i.i.d. Gaussian. 



Robust PCA. PCA is arguably the most widely used statistical tool for data analysis and dimen- 
sionality reduction today. However, its brittleness with respect to grossly corrupted observations 
often puts its validity in jeopardy - a single grossly corrupted entry in M could render the estimated 
L arbitrarily far from the true Lq. Unfortunately, gross errors are now ubiquitous in modern appli- 
cations such as image processing, web data analysis, and bioinformatics, where some measurements 
may be arbitrarily corrupted (due to occlusions, malicious tampering, or sensor failures) or simply 
irrelevant to the low-dimensional structure we seek to identify. A number of natural approaches 
to robustifying PCA have been explored and proposed in the literature over several decades. The 
representative approaches include influence function techniques [26,47], multivariate trimming [19], 
alternating minimization [28], and random sampling techniques [17]. Unfortunately, none of these 
existing approaches yields a polynomial-time algorithm with strong performance guarantees under 
broad conditions^] The new problem we consider here can be considered as an idealized version of 
Robust PCA, in which we aim to recover a low-rank matrix Lq from highly corrupted measurements 
M = Lq + Sq. Unlike the small noise term Nq in classical PCA, the entries in Sq can have arbitrarily 
large magnitude, and their support is assumed to be sparse but unknown^} 

2 We refer to either the complexity of algorithms that increases drastically as dimension increases, or to their 
performance that decreases sharply when scale goes up. 

3 Random sampling approaches guarantee near-optimal estimates, but have complexity exponential in the rank of 
the matrix Lo. Trimming algorithms have comparatively lower computational complexity, but guarantee only locally 
optimal solutions. 

4 The unknown support of the errors makes the problem more difficult than the matrix completion problem that 
has been recently much studied. 
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Applications. There are many important applications in which the data under study can natu- 
rally be modeled as a low-rank plus a sparse contribution. All the statistical applications, in which 
robust principal components are sought, of course fit our model. Below, we give examples inspired 
by contemporary challenges in computer science, and note that depending on the applications, 
either the low-rank component or the sparse component could be the object of interest: 

• Video Surveillance. Given a sequence of surveillance video frames, we often need to identify 
activities that stand out from the background. If we stack the video frames as columns 
of a matrix M, then the low-rank component L$ naturally corresponds to the stationary 
background and the sparse component Sq captures the moving objects in the foreground. 
However, each image frame has thousands or tens of thousands of pixels, and each video 
fragment contains hundreds or thousands of frames. It would be impossible to decompose M 
in such a way unless we have a truly scalable solution to this problem. In Section |4j we will 
show the results of our algorithm on video decomposition. 

• Face Recognition. It is well known that images of a convex, Lambertian surface under varying 
illuminations span a low-dimensional subspace [1]. This fact has been a main reason why 
low-dimensional models are mostly effective for imagery data. In particular, images of a 
human's face can be well-approximated by a low-dimensional subspace. Being able to correctly 
retrieve this subspace is crucial in many applications such as face recognition and alignment. 
However, realistic face images often suffer from self-shadowing, specularities, or saturations 
in brightness, which make this a difficult task and subsequently compromise the recognition 
performance. In Section |4j we will show how our method is able to effectively remove such 
defects in face images. 

• Latent Semantic Indexing. Web search engines often need to analyze and index the content 
of an enormous corpus of documents. A popular scheme is the Latent Semantic Indexing 
(LSI) [14,42]. The basic idea is to gather a document-versus-term matrix M whose entries 
typically encode the relevance of a term (or a word) to a document such as the frequency it 
appears in the document (e.g. the TF/IDF). PCA (or SVD) has traditionally been used to 
decompose the matrix as a low-rank part plus a residual, which is not necessarily sparse (as 
we would like) . If we were able to decompose M as a sum of a low- rank component Lq and a 
sparse component So, then Lq could capture common words used in all the documents while 
Sq captures the few key words that best distinguish each document from others. 

• Ranking and Collaborative Filtering. The problem of anticipating user tastes is gaining in- 
creasing importance in online commerce and advertisement. Companies now routinely collect 
user rankings for various products, e.g., movies, books, games, or web tools, among which 
the Netflix Prize for movie ranking is the best known [40] . The problem is to use incomplete 
rankings provided by the users on some of the products to predict the preference of any given 
user on any of the products. This problem is typically cast as a low-rank matrix completion 
problem. However, as the data collection process often lacks control or is sometimes even 
ad hoc - a small portion of the available rankings could be noisy and even tampered with. 
The problem is more challenging since we need to simultaneously complete the matrix and 
correct the errors. That is, we need to infer a low-rank matrix Lq from a set of incomplete 



and corrupted entries. In Section 1.6 we will see how our results can be extended to this 
situation. 
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Similar problems also arise in many other applications such as graphical model learning, linear 
system identification, and coherence decomposition in optical systems, as discussed in [12]. All in 
all, the new applications we have listed above require solving the low-rank and sparse decomposition 
problem for matrices of extremely high dimension and under much broader conditions, a goal this 
paper aims to achieve. 



1.2 A surprising message 

At first sight, the separation problem seems impossible to solve since the number of unknowns to 
infer for Lq and So is twice as many as the given measurements in M E M niXn2 . Furthermore, it 
seems even more daunting that we expect to reliably obtain the low-rank matrix Lq with errors in 
So of arbitrarily large magnitude. 

In this paper, we are going to see that very surprisingly, not only can this problem be solved, 
it can be solved by tractable convex optimization. Let := Y2i a i(M) denote the nuclear 

norm of the matrix M, i.e. the sum of the singular values of M, and let ||M||i = £^ |-&%| denote 
the £i-norm of M seen as a long vector in IR niXTl2 . Then we will show that under rather weak 
assumptions, the Principal Component Pursuit (PCP) estimate solving^ 



minimize + A||S||i , , 

subject to L + S = M ^ ' ' 

exactly recovers the low-rank Lq and the sparse So- Theoretically, this is guaranteed to work even if 
the rank of Lq grows almost linearly in the dimension of the matrix, and the errors in Sq are up to a 
constant fraction of all entries. Algorithmically, we will see that the above problem can be solved by 
efficient and scalable algorithms, at a cost not so much higher than the classical PCA. Empirically, 
our simulations and experiments suggest this works under surprisingly broad conditions for many 



types of real data. In Section 1.5 we will comment on the similar approach taken in the paper [12], 



which was released during the preparation of this manuscript. 



1.3 When does separation make sense? 

A normal reaction is that the objectives of this paper cannot be met. Indeed, there seems to not be 
enough information to perfectly disentangle the low-rank and the sparse components. And indeed, 
there is some truth to this, since there obviously is an identifiability issue. For instance, suppose 
the matrix M is equal to e\e\ (this matrix has a one in the top left corner and zeros everywhere 
else). Then since M is both sparse and low-rank, how can we decide whether it is low-rank or 
sparse? To make the problem meaningful, we need to impose that the low-rank component Lq is 
not sparse. In this paper, we will borrow the general notion of incoherence introduced in [8] for the 
matrix completion problem; this is an assumption concerning the singular vectors of the low-rank 
component. Write the singular value decomposition of Lq £ M n i x ™2 as 

r 

Lq = UEV* = ^ OiUiV* , 
i=l 

5 Although the name naturally suggests an emphasis on the recovery of the low-rank component, we reiterate that 
in some applications, the sparse component truly is the object of interest. 
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where r is the rank of the matrix, a\, . . . ,o~ r are the positive singular values, and U = [ui, . . . , u r ], 
V = [vi , . . . , v r ] are the matrices of left- and right-singular vectors. Then the incoherence condition 
with parameter p states that 

Wtt* ||2 ^ ^ r NT/* \\2 ^ / ir /I n \ 

max t; e$ < — , max k ed < — , (1-2) 
i mi n ■> 

and 



Halloo < J — • (1.3) 

V n l n 2 

Here and below, HMH^ = maxjj |M«|, i.e. is the £oo norm of M seen as a long vector. Note 
that since the orthogonal projection Pjj onto the column space of U is given by Pjj = UU* , (1.2) 
is equivalent to max« ||P(/ej|| 2 < pr/ni, and similarly for Py. As discussed in earlier references 
[8,10,22], the incoherence condition asserts that for small values of p, the singular vectors are 
reasonably spread out - in other words, not sparse. 

Another identifiability issue arises if the sparse matrix has low-rank. This will occur if, say, all 
the nonzero entries of S occur in a column or in a few columns. Suppose for instance, that the first 
column of So is the opposite of that of Lq, and that all the other columns of So vanish. Then it is 
clear that we would not be able to recover Lq and So by any method whatsoever since M = Lq + So 
would have a column space equal to, or included in that of Lq. To avoid such meaningless situations, 
we will assume that the sparsity pattern of the sparse component is selected uniformly at random. 

1.4 Main result 

The surprise is that under these minimal assumptions, the simple PCP solution perfectly recovers 
the low- rank and the sparse components, provided of course that the rank of the low- rank compo- 
nent is not too large, and that the sparse component is reasonably sparse. Below, nt\\ = max(ni, 77-2) 
and n(2) = min(ni, 712). 



Theorem 1.1 Suppose Lq is n x n, obeys (1.2)-(1.3), and that the support set of Sq is uniformly 



distributed among all sets of cardinality m. Then there is a numerical constant c such that with 
probability at least 1 — cn -10 ( over the choice of support of Sq ), Principal Component Pursuit ( 1.1 ) 
with A = 1/y/n is exact, i.e. L = Lq and S = Sq, provided that 

rank(Lo) < p r n /i _1 (log n) -2 and m<p s n 2 . (1-4) 

Above, p r and p s are positive numerical constants. In the general rectangular case where Lq is 
n\XU2, PCP with A = 1/ succeeds with probability at least 1 — cnZy-P, provided that rank(Lo) < 

/5 r n( 2 ) A i_1 (l°g n (i)) _2 and m < p s n\U2- 

In other words, matrices Lq whose singular vectors — or principal components — are reasonably 
spread can be recovered with probability nearly one from arbitrary and completely unknown cor- 
ruption patterns (as long as these are randomly distributed). In fact, this works for large values of 
the rank, i.e. on the order of n/(logn) 2 when p is not too large. We would like to emphasize that 
the only 'piece of randomness' in our assumptions concerns the locations of the nonzero entries 
of So; everything else is deterministic. In particular, all we require about Lq is that its singular 
vectors are not spiky. Also, we make no assumption about the magnitudes or signs of the nonzero 
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entries of So- To avoid any ambiguity, our model for So is this: take an arbitrary matrix S and set 
to zero its entries on the random set O c ; this gives Sq. 

A rather remarkable fact is that there is no tuning parameter in our algorithm. Under the 
assumption of the theorem, minimizing 

\\L\\* H ll^lli; n (D = maxfrii, n 2 ) 

always returns the correct answer. This is surprising because one might have expected that one 
would have to choose the right scalar A to balance the two terms in ||L||* + A||5||i appropriately 
(perhaps depending on their relative size). This is, however, clearly not the case. In this sense, the 
choice A = 1/ is universal. Further, it is not a priori very clear why A = 1/ ^/n(i) is a correct 
choice no matter what Lq and Sq are. It is the mathematical analysis which reveals the correctness 
of this value. In fact, the proof of the theorem gives a whole range of correct values, and we have 
selected a sufficiently simple value in that range. 

Another comment is that one can obtain results with larger probabilities of success, i.e. of the 
form 1 — 0(n~P) (or 1 — 0(n^)) for [3 > at the expense of reducing the value of p r . 

1.5 Connections with prior work and innovations 

The last year or two have seen the rapid development of a scientific literature concerned with the 
matrix completion problem introduced in [8], see also [7,10,22,23,43] and the references therein. In 
a nutshell, the matrix completion problem is that of recovering a low- rank matrix from only a small 
fraction of its entries, and by extension, from a small number of linear functionals. Although other 
methods have been proposed [43], the method of choice is to use convex optimization [7,10,22,23,45]: 
among all the matrices consistent with the data, simply find that with minimum nuclear norm. 
The papers cited above all prove the mathematical validity of this approach, and our mathematical 
analysis borrows ideas from this literature, and especially from those pioneered in [8] . Our methods 
also much rely on the powerful ideas and elegant techniques introduced by David Gross in the 
context of quantum-state tomography [22,23]. In particular, the clever golfing scheme [22] plays a 
crucial role in our analysis, and we introduce two novel modifications to this scheme. 

Despite these similarities, our ideas depart from the literature on matrix completion on several 
fronts. First, our results obviously are of a different nature. Second, we could think of our sep- 
aration problem, and the recovery of the low-rank component, as a matrix completion problem. 
Indeed, instead of having a fraction of observed entries available and the other missing, we have 
a fraction available, but do not know which one, while the other is not missing but entirely cor- 
rupted altogether. Although, this is a harder problem, one way to think of our algorithm is that 
it simultaneously detects the corrupted entries, and perfectly fits the low-rank component to the 
remaining entries that are deemed reliable. In this sense, our methodology and results go beyond 
matrix completion. Third, we introduce a novel de-randomization argument that allows us to fix 
the signs of the nonzero entries of the sparse component. We believe that this technique will have 
many applications. One such application is in the area of compressive sensing, where assumptions 
about the randomness of the signs of a signal are common, and merely made out of convenience 
rather than necessity; this is important because assuming independent signal signs may not make 
much sense for many practical applications when the involved signals can all be non-negative (such 
as images). 
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We mentioned earlier the related work [12], which also considers the problem of decomposing a 
given data matrix into sparse and low-rank components, and gives sufficient conditions for convex 
programming to succeed. These conditions are phrased in terms of two quantities. The first is the 
maximum ratio between the ^oo norm and the operator norm, restricted to the subspace generated 
by matrices whose row or column spaces agree with those of Lq. The second is the maximum ratio 
between the operator norm and the norm, restricted to the subspace of matrices that vanish 
off the support of Sq. Chandrasekaran et. al. show that when the product of these two quantities 
is small, then the recovery is exact for a certain interval of the regularization parameter [12]. 

One very appealing aspect of this condition is that it is completely deterministic: it does not 
depend on any random model for Lq or So- It yields a corollary that can be easily compared to 
our result: suppose n\ = n<i = n for simplicity, and let hq be the smallest quantity satisfying (1.2), 
then correct recovery occurs whenever 

max{i : [So]ij / 0} x \J ^xtfrjn < 1/12. 

3 

The left-hand side is at least as large as p sy /jMynr, where p s is the fraction of entries of So that are 
nonzero. Since fiQ > 1 always, this statement only guarantees recovery if p s = 0((nr)~ 1 / 2 ); i.e., 
even when rank(Lo) = 0(1), only vanishing fractions of the entries in Sq can be nonzero. 

In contrast, our result shows that for incoherent Lq, correct recovery occurs with high probability 
for rank(Lo) on the order of n/[/ilog 2 n] and a number of nonzero entries in Sq on the order of n 2 . 
That is, matrices of large rank can be recovered from non- vanishing fractions of sparse errors. This 
improvement comes at the expense of introducing one piece of randomness: a uniform model on 
the error support 

Our analysis has one additional advantage, which is of significant practical importance: it iden- 
tifies a simple, non-adaptive choice of the regularization parameter A. In contrast, the conditions 
on the regularization parameter given by Chandrasekaran et al. depend on quantities which in 
practice are not known a-priori. The experimental section of [12] suggests searching for the correct 
A by solving many convex programs. Our result, on the other hand, demonstrates that the simple 
choice A = 1 j \fn works with high probability for recovering any square incoherent matrix. 

1.6 Implications for matrix completion from grossly corrupted data 

We have seen that our main result asserts that it is possible to recover a low-rank matrix even 
though a significant fraction of its entries are corrupted. In some applications, however, some 
of the entries may be missing as well, and this section addresses this situation. Let Vq be the 
orthogonal projection onto the linear space of matrices supported on f2 C [ni] x [ri2\ , 



Then imagine we only have available a few entries of Lq + Sq, which we conveniently write as 

Y = Va obB (Lo + So)=Vn obB LQ + S'o; 



6 Notice that the bound of [12] depends only on the support of So, and hence can be interpreted as a worst case 
result with respect to the signs of So- In contrast, our result does not randomize over the signs, but does assume that 
they are sampled from a fixed sign pattern. Although we do not pursue it here due to space limitations, our analysis 
also yields a result which holds for worst case sign patterns, and guarantees correct recovery with rank(Lo) = 0(1), 
and a sparsity pattern of cardinality pn\ni for some p > 0. 
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that is, we see only those entries G S^obs C [ni] x [742]- This models the following problem: we 
wish to recover Lq but only see a few entries about Lq, and among those a fraction happens to be 
corrupted, and we of course do not know which one. As is easily seen, this is a significant extension 
of the matrix completion problem, which seeks to recover Lq from undersampled but otherwise 
perfect data T>n ohs L . 

We propose recovering Lq by solving the following problem: 

minimize HLH* + A||5||i , , 

subject to Vn ohs (L + S)=Y. ( - 5j 

In words, among all decompositions matching the available data, Principal Component Pursuit 
finds the one that minimizes the weighted combination of the nuclear norm, and of the l\ norm. Our 
observation is that under some conditions, this simple approach recovers the low- rank component 
exactly. In fact, the techniques developed in this paper establish this result: 



Theorem 1.2 Suppose Lq is n x n, obeys the conditions (1.2)-(1.3), and that &s is uniformly 
distributed among all sets of cardinality m obeying m = O.Ira 2 . Suppose for simplicity, that each 
observed entry is corrupted with probability r independently of the others. Then there is a numerical 
constant c such that with probability at least 1 — cn -10 , Principal Component Pursuit (1.5 1 with 
A = 1/yO.ln is exact, i.e. L = Lq, provided that 

rank(Lo) < p r rau _1 (logra) -2 , and r < r s . (1-6) 

Above, p r and t s are positive numerical constants. For general m x n 2 rectangular matrices, PCP 
with A = l/^O.lra^) succeeds from m = O.lnira.2 corrupted entries with probability at least 1 — cra^°, 

provided that rank(Lo) < p 7 .7i(2)A' (lognm) -2 . 

In short, perfect recovery from incomplete and corrupted entries is possible by convex optimization. 

On the one hand, this result extends our previous result in the following way. If all the entries 
are available, i.e. m = n\n2, then this is Theorem On the other hand, it extends matrix 
completion results. Indeed, if r = 0, we have a pure matrix completion problem from about a 
fraction of the total number of entries, and our theorem guarantees perfect recovery as long as r 



obeys (1.6), which for large values of r, matches the strongest results available. We remark that 
the recovery is exact, however, via a different algorithm. To be sure, in matrix completion one 
typically minimizes the nuclear norm subject to the constraint Vn obs L = Vn oha LQ. Here, our 
program would solve 

minimize \\L\\* + A||5||i (17) 
subject to Vn obs (L + S) = Vn obs L , ( ' ' 



and return L = Lq, S = 0! In this context, Theorem 1.2 proves that matrix completion is stable 
vis a vis gross errors. 



Remark. We have stated Theorem 1.2 merely to explain how our ideas can easily be adapted 
to deal with low-rank matrix recovery problems from undersampled and possibly grossly corrupted 
data. In our statement, we have chosen to see 10% of the entries but, naturally, similar results 
hold for all other positive fractions provided that they are large enough. We would like to make it 



clear that a more careful study is likely to lead to a stronger version of Theorem 1.2 In particular 



for very low rank matrices, we expect to see similar results holding with far fewer observations; 
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that is, in the limit of large matrices, from a decreasing fraction of entries. In fact, our techniques 
would already establish such sharper results but we prefer not to dwell on such refinements at the 
moment, and leave this up for future work. 



1.7 Notation 

We provide a brief summary of the notations used throughout the paper. We shall use five norms 
of a matrix. The first three are functions of the singular values and they are: 1) the operator norm 
or 2-norm denoted by \\X\\; 2) the Frobenius norm denoted by and 3) the nuclear norm 

denoted by ||^||*. The last two are the l\ and 1^ norms of a matrix seen as a long vector, and 
are denoted by ||-X"||i and ||^||oo respectively. The Euclidean inner product between two matrices 
is defined by the formula (X, Y) := trace(X*Y), so that ||X|||, = (X,X). 

Further, we will also manipulate linear transformations which act on the space of matrices, and 
we will use calligraphic letters for these operators as in Vq.X. We shall also abuse notation by also 
letting £1 be the linear space of matrices supported on 0,. Then T > q± denotes the projection onto 
the space of matrices supported on £l c so that X = Vq + Vq± , where X is the identity operator. We 
will consider a single norm for these, namely, the operator norm (the top singular value) denoted 
by ||^4||, which we may want to think of as ||^4|| = sup/n^ii^n ||.4X||_f; for instance, \\Vq\\ = 1 
whenever 0, ^ 0. 



1.8 Organization of the paper 

The paper is organized as follows. In Section [2] we provide the key steps in the proof of Theorem 
This proof depends upon on two critical properties of dual certificates, which are established in 
the separate Section [3} The reason why this is separate is that in a first reading, the reader might 
want to jump to Section |4j which presents applications to video surveillance, and computer vision. 
Section [5] introduces algorithmic ideas to find the Principal Component Pursuit solution when M 
is of very large scale. We conclude the paper with a discussion about future research directions in 
Section [6] Finally, the proof of Theorem 1.2 is in the Appendix, Section [7j together with those of 
intermediate results. 



2 Architecture of the Proof 

This section introduces the key steps underlying the proof of our main result, Theorem We 
will prove the result for square matrices for simplicity, and write n = n\ = rii- Of course, we shall 
indicate where the argument needs to be modified to handle the general case. Before we start, 
it is helpful to review some basic concepts and introduce additional notation that shall be used 
throughout. For a given scalar x, we denote by sgn(x) the sign of x, which we take to be zero if 
x = 0. By extension, sgn(5) is the matrix whose entries are the signs of those of S. We recall that 
any subgradient of the l\ norm at Sq supported on Q, is of the form 

sgn(,So) + F, 

where F vanishes on f2, i.e. VqF = 0, and obeys ||-F||oo ^ 1- 

We will also manipulate the set of subgradients of the nuclear norm. From now on, we will 
assume that Lq of rank r has the singular value decomposition UT,V*, where U, V € W nxr just as 
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in Section 1.3 Then any subgradient of the nuclear norm at Lq is of the form 

uv* + w, 

where U*W = 0, WV = and ||W|| < 1. Denote by T the linear space of matrices 

T := {UX* + YV*, X, Y £ R nxr }, (2.1) 

and by T 1 - its orthogonal complement. It is not hard to see that taken together, U*W = and 
WV = are equivalent to VtW = 0, where Vt is the orthogonal projection onto T. Another way to 
put this is V T ±W = W. In passing, note that for any matrix M, V T± M = (I - UU*)M(I - VV*), 
where we recognize that / — UU* is the projection onto the orthogonal complement of the linear 
space spanned by the columns of U and likewise for (I — VV*). A consequence of this simple 
observation is that for any matrix M, \\Vt±M\\ < \\M\\, a fact that we will use several times in the 
sequel. Another consequence is that for any matrix of the form e^e*, 

\\V T ±e ie *\\ 2 F = ||(/ - UU*)ei\\ 2 \\(I - VV*)e 3 f > (1 - fir/n) 2 , 

where we have assumed jir/n < 1. Since WPx^i&j \\ 2 f + WPT ±e i e *j [If = 1) t ms gives 

\\V T e l e)\\ F <J 2 ^. (2.2) 



For rectangular matrices, the estimate is [["PT^e* ||f < J m i n (^ n2 ) • 

Finally, in the sequel we will write that an event holds with high or large probability whenever 

Hi) 



it holds with probability at least 1 — 0(n 10 ) (with nn\ in place of n for rectangular matrices). 



2.1 An elimination theorem 

We begin with a useful definition and an elementary result we shall use a few times. 

Definition 2.1 We will say that S' is a trimmed version of S if supp(S') C supp(S) and S'^ = Sij 
whenever S'^ / 0. 

In words, a trimmed version of S is obtained by setting some of the entries of S to zero. Having said 
this, the following intuitive theorem asserts that if Principal Component Pursuit correctly recovers 
the low-rank and sparse components of Mq = Lq + So, it also correctly recovers the components of 
a matrix M = Lq + S' where S' Q is a trimmed version of So- This is intuitive since the problem is 
somehow easier as there are fewer things to recover. 



Theorem 2.2 Suppose the solution to (1.1) with input data Mq = Lq + So is unique and exact, 
and consider Mq = Lq + S' , where S' is a trimmed version of So. Then the solution to (1.1 ) with 
input Mq is exact as well. 



Proof Write S' = V^Sq for some Qo C [n] x [n] and let (L, S) be the solution of ( 1.1 ) with input 
Lq + S' . Then 

||L||* + A||5||i < ||Lo||* + A||7>n So[|i 

and, therefore, 

||£||* + A||5||i + A||P n ±5 ||i < \\L \U + X\\So\\i. 
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Note that (L,S + VqxSo) is feasible for the problem with input data Lq + So, and since \\S + 
TV 'So 111 ^ ll*S||i + ||^o-i5 ||i, we have 



\\L\U + \\\S + V n ±S \\i < ||L ||* + A||5o||i. 

The right-hand side, however, is the optimal value, and by unicity of the optimal solution, we must 
have L = Lq, and S + V n ±So = So or S = Vq So = S' Q . This proves the claim. ■ 

The Bernoulli model. In Theorem probability is taken with respect to the uniformly 
random subset 0, = ■ S{j ^ 0} of cardinality m. In practice, it is a little more convenient to 

work with the Bernoulli model £1 = '■ Sij = 1}, where the dij's are i.i.d. variables Bernoulli 

taking value one with probability p and zero with probability 1 — p, so that the expected cardinality 
of $7 is pn 2 . From now on, we will write f2 ~ Ber(/?) as a shorthand for f2 is sampled from the 
Bernoulli model with parameter p. 



Since by Theorem 2.2 the success of the algorithm is monotone in |0| , any guarantee proved for 
the Bernoulli model holds for the uniform model as well, and vice versa, if we allow for a vanishing 
shift in p around m/n 2 . The arguments underlying this equivalence are standard, see [9, 10], and 
may be found in the Appendix for completeness. 

2.2 Derandomization 

In Theorem the values of the nonzero entries of So are fixed. It turns out that it is easier 
to prove the theorem under a stronger assumption, which assumes that the signs of the nonzero 
entries are independent symmetric Bernoulli variables, i.e. take the value ±1 with probability 
1/2 (independently of the choice of the support set). The convenient theorem below shows that 
establishing the result for random signs is sufficient to claim a similar result for fixed signs. 



Theorem 2.3 Suppose Lq obeys the conditions of Theorem 1.1 and that the locations of the nonzero 
entries of Sq follow the Bernoulli model with parameter 2p s , and the signs of So are i.i.d. ±1 as 
above (and independent from the locations). Then if the PCP solution is exact with high probability, 
then it is also exact with at least the same probability for the model in which the signs are fixed and 
the locations are sampled from the Bernoulli model with parameter p s . 

This theorem is convenient because to prove our main result, we only need to show that it is true 
in the case where the signs of the sparse component are random. 

Proof Consider the model in which the signs are fixed. In this model, it is convenient to think of 
So as VqS, for some fixed matrix S, where f2 is sampled from the Bernoulli model with parameter 
p s . Therefore, So has independent components distributed as 



(So, 



ij 



Sij, W. p. p s , 

0, w. p. 1 - p s . 



Consider now a random sign matrix with i.i.d. entries distributed as 



Eij 



1, w. p. p a , 

0, w. p. 1 — 2p s 

^-1, W. p. p s , 
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and an "elimination" matrix A with entries denned by 

A fo, if^[sgn(S)]y = -1, 
[l, otherwise. 

Note that the entries of A are independent since they are functions of independent variables. 

Consider now S' Q = A o (|S| o E), where o denotes the Hadamard or componentwise product so 
that, [S'o\ij = Ajj (\Sij\Eij). Then we claim that S' Q and So have the same distribution. To see why 
this is true, it suffices by independence to check that the marginals match. For Sij ^ 0, we have 

P([S&]y = Sij) = P(A ii = 1 and Eij = [sgn(S)y 

= P(^[sgn(5)]y + -1 and E i5 = [sgn(S)]y) 
= F(Eij = [sgn(S)k) = Ps, 

which establishes the claim. 

This construction allows to prove the theorem. Indeed, |S|oL now obeys the random sign model, 
and by assumption, PCP recovers |S| o E with high probability. By the elimination theorem, this 
program also recovers S' Q = A o (|S| o E). Since S' and So have the same distribution, the theorem 
follows. ■ 



2.3 Dual certificates 

We introduce a simple condition for the pair (Lo, So) to be the unique optimal solution to Principal 
Component Pursuit. These conditions are stated in terms of a dual vector, the existence of which 
certifies optimality. (Recall that O is the space of matrices with the same support as the sparse 
component So, and that T is the space defined via the the column and row spaces of the low-rank 
component Lo (2.1).) 

Lemma 2.4 Assume that \\VqVt\\ < 1> With the standard notations, (Lo,Sq) is the unique solu- 
tion if there is a pair (W, F) obeying 

UV* + W = X(sgn(S ) + F), 

with V T W = 0, ||W|| < I, V n F = and WF]]^ < 1. 

Note that the condition H'Pq'PtII < 1 is equivalent to saying that O H T = {0}. 
Proof We consider a feasible perturbation (Lo + H, So — H) and show that the objective increases 
whenever H ^ 0, hence proving that (Lo, So) is the unique solution. To do this, let UV* + Wo be 
an arbitrary subgradient of the nuclear norm at Lo, and sgn(So) + Lb be an arbitrary subgradient 
of the ^i-norm at So- By definition of subgradients, 

||L + H\U + A|| So - # 111 > ll^oll* + A||S ||i + (UV* + W , H) - A(sgn(S ) + L , H). 

Now pick Wo such that (W ,H) = HTVlL^I* and L such that {F ,H) = -\\P n ±H\\i^\ We have 

||L + F||, + A||S - H\\i > \\L \U + A||5 ||i + ll^^ll* + Ml^H^ + (UV* - Asgn(S ), H). 



7 For instance, Fq = — sgn(V n ± H) is such a matrix. Also, by duality between the nuclear and the operator norm, 
there is a matrix obeying || W|| = 1 such that (W,V T ±H) — \\V T ±H\\*, and we just take Wo = V T ± (W). 
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By assumption 

\(UV* -\ S gn(S ),H}\ < \(W,H)\ + X\(F,H)\ < + APV#lli) 
for (3 = max(||W||, Halloo) < 1 and, thus, 

\\Lq + H\U + X\\S - > ||f-o||* + A||5 ||i + (1-0) (\\v t ±h\u + a||tv#IIi 

Since by assumption, f2 n T = {0}, we have || Pyx fry* + AHP^xifHi > unless H = 0. ■ 

Hence, we see that to prove exact recovery, it is sufficient to produce a 'dual certificate' W 
obeying 

'w e T L , 
W\\ < l, 

P Q ([/y* + W0 = Asgn(So), l2 " >) 

j|p n x(rjy* + w)||oo<A. 

Our method, however, will produce with high probability a slightly different certificate. The idea 
is to slightly relax the constraint Pq,(UV* + W) = Asgn(So), a relaxation that has been introduced 
by David Gross in [22] in a different context. We prove the following lemma. 

Lemma 2.5 Assume WVqPtW < 1/2 and A < 1. Then with the same notation, (Lq,Sq) is the 
unique solution if there is a pair (W, F) obeying 

UV* + W = X(sgn(S ) + F + P a D) 

with P T W = and \\W\\ < \, P n F = and WF]]^ < \, and \\PqD\\ f < \. 



Proof Following the proof of Lemma 2.4 we have 

||f o + H\\, + X\\S - fT||i > ||Lo||* + A||5 ||i + i (\\V T ±H\\* + X\\V n ±H\\^j - X{P n D, H) 

> \\L \U + A||5 ||i + \ (||7V# II* + A||7V# Hi) - ^ll^nfTllF- 



Observe now that 



and, therefore, 
In conclusion, 



\\PqH\\ f < \\V n V T H\\ F + \\PnP T xH\\ F 
— 2 ll-^ll-F \\T~ > t ± H\\f 
< l -\\VnH\\ F + ^\\V Q xH\\ F + \\V T ±H\\ F 



\\VnH\\ F < \\P a ±H\\ F + 2\\P T xH\\ F . 



\Lq + h\u + x\\s - fT||i > ||f oil* + A||5 ||i + l((i- X)\[P T xH\U + ^||7Vfr||i' 
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and the term between parenthesis is strictly positive when H ^ 0. 

As a consequence of Lemma |2.5| it now suffices to produce a dual certificate W obeying 



W £ T , 

\\w\\ < 1/2, 

\\V n (UV* - Asgn(5 ) + W)\\ F < A/4, 
J|7V([/F* + W0||oo<A/2. 

Further, we would like to note that the existing literature on matrix completion [8] gives good 



bounds on ||7 ? q7 ? t|| j see Theorem 2.6 in Section 2.5 



2.4 Dual certification via the golfing scheme 

In the papers [22,23], Gross introduces a new scheme, termed the golfing scheme, to construct a 
dual certificate for the matrix completion problem, i.e. the problem of reconstructing a low-rank 
matrix from a subset of its entries. In this section, we will adapt this clever golfing scheme, with 
two important modifications, to our separation problem. 

Before we introduce our construction, our model assumes that f2 ~ Ber(p), or equivalently that 
f2 c ~ Ber(l — p). Now the distribution of £l c is the same as that of £l c = 0,± U ^2 U . . . U fL- , where 
each £lj follows the Bernoulli model with parameter q, which has an explicit expression. To see 
this, observe that by independence, we just need to make sure that any entry is selected with 
the right probability. We have 

P((i, 3) £«) = P(Bin(j ,<z) = 0) = (1 - qy\ 
so that the two models are the same if 

P = (i- q y°, 

hence justifying our assertion. Note that because of overlaps between the Oj's, q > (1 — p)/jo- 
We now propose constructing a dual certificate 

W = W L + W s , 

where each component is as follows: 

1. Construction of W L via the golfing scheme. Fix an integer jo > 1 whose value shall be 
discussed later, and let Qj, 1 < j < jo, be defined as above so that fl c = ^->i<j<j ^j- Then 
starting with Yq = 0, inductively define 

and set 

W L = V T± Y j0 . (2.5) 

This is a variation on the golfing scheme discussed in [22], which assumes that the fJj's 
are sampled with replacement, and does not use the projector *Pq. but something more 
complicated taking into account the number of times a specific entry has been sampled. 
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2. Construction of W s via the method of least squares. Assume that (('P^'PtII < 1/2. Then 
WPn'Pr'PnW < 1/4 and, thus, the operator Vn — VqVtVq mapping f2 onto itself is invertible; 
we denote its inverse by (Pa - VuVrVny 1 . We then set 

W s = XV T ±(Vn - VnVTVnr^gniSo). (2.6) 

Clearly, an equivalent definition is via the convergent Neumann series 

W s = AP T x^(PnPT^) fc sgn( 1 S ). (2.7) 

k>0 

Note that V Q W S = XP n (I - V T )(V n - VoPtPq)- 1 ^^) = Asgn(S ). With this, the 
construction has a natural interpretation: one can verify that among all matrices W £ T 1 - 
obeying VqW = Asgn(So), W s is that with minimum Frobenius norm. 

Since both W L and W s belong to T 1 - and VqW s = Asgn(5 ), we will establish that W L + W s is 
a valid dual certificate if it obeys 

\\W L + W s \\ < 1/2, 

\\V n (UV* + W L )\\ F <\/4, (2.8) 

[\\r n ±(uv* + w L + w s )\\ O0 < A/2. 

2.5 Key lemmas 

We now state three lemmas, which taken collectively, establish our main theorem. The first may 
be found in [8] . 

Theorem 2.6 [8, Theorem 4-1] Suppose is sampled from the Bernoulli model with parameter 
Pq. Then with high probability, 

\\V T ~ P^VtV^VtW <e, (2.9) 

provided that p > Co e" 2 for some numerical constant Co > (p is the incoherence param- 

eter). For rectangular matrices, we need po > Co e~ 2 • 

Among other things, this lemma is important because it shows that HPhPtII < 1/2, provided 
|0| is not too large. Indeed, if ~ Ber(p), we have 

\\V T - (1 - pY^tV^VtW <e, 
with the proviso that 1 — p > Co e~ 2 Mr '° g - . Note, however, that since 1 = Vn + V n ±, 

V T -(1- pT^tV^Vt = (1 - pyHV T VnVT - pV T ) 
and, therefore, by the triangular inequality 

WVtVqVtW < e(l -p) + p\\V T \\ =p + e(l- p). 
Since HTVPtII 2 = \\PtVqPt \\- we have established the following: 
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Corollary 2.7 Assume that f2 ~ Ber(p), then \\VqVt\\ 2 < p+e, provided that 1—p > Cq e 



-2 Uri 



where Cq is as in Theorem 2.6. For rectangular matrices, the modification is as in Theorem 2.6 



The lemma below is proved is Section [3] 

Lemma 2.8 Assume that VL ~ Ber{p) with parameter p < p s for some p s > 0. Set jo = 2[~logn] 
(use logn(x) for rectangular matrices). Then under the other assumptions of Theorem the 
matrix W L (|2.5|) obeys 



(a) \\W L \\ < 1/4, 

(b) \\Vn{UV* + W L )\\ F < A/4, 

(c) \\V^{UV* + W L )\\ 00 <X/A. 

Since H'Pfi'PTll < 1 with large probability, W s is well defined and the following holds. 



Lemma 2.9 Assume that Sq is supported on a set $7 sampled as in Lemma 2.8, and that the signs 
of Sq are i.i.d. symmetric (and independent of£l). Then under the other assumptions of Theorem 
\l.l\ the matrix W s ( |2.6[ ) obeys 

(a) \\W S \\ < 1/4, 

(b) ||7VW 5 ||oo < V 4 - 

The proof is also in Section 3 Clearly, W L and W obey (2.8), hence certifying that Principal 



Component Pursuit correctly recovers the low-rank and sparse components with high probabil- 
ity when the signs of Sq are random. The earlier "derandomization" argument then establishes 
Theorem 11.11 



3 Proofs of Dual Certification 



This section proves the two crucial estimates, namely, Lemma 2.8 and Lemma 2.9 



3.1 Preliminaries 



We begin by recording two results which shall be useful in proving Lemma 2.8 While Theorem 2.6 
asserts that with large probability, 

\\ z ~ Pq 1 'Pt'Pq, Z\\ f < e\\Z\\ F , 

for all Z G T, the next lemma shows that for a fixed Z, the sup-norm of Z — p ~ 1 Vt'Pq (Z) also 
does not increase (also with large probability). 

Lemma 3.1 Suppose Z £ T is a fixed matrix, and 0,q ~ Ber{po). Then with high probability, 

\\Z - ps l V T Vn Q Z\\oo < eWZW^ (3.1) 
provided that po > Cq e~ 2 ^ r log - (for rectangular matrices, po > Co e -2 ^ r lo g n a) j j or some numer- 

(2) 

ical constant Cq > 0. 
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The proof is an application of Bernstein's inequality and may be found in the Appendix. A similar 
but somewhat different version of (3.1) appears in [44]. 
The second result was proved in [8]. 

Lemma 3.2 [8, Theorem 6.3] Suppose Z is fixed, and £Iq ~ Ber(po). Then with high probability, 

IKJ-Po 1 ^)^!! <Co 



' n log n , 



Z 



PO 



(3.2) 



for some small numerical constant C' Q > provided that po > Co - ^ g "' (or po > C' Atl ° gn ( 1 ) j or 
rectangular matrices in which case nm log nm replaces relogra in (3.2)). 



As a remark, Lemmas 



3.1 



and 



3.2 



and Theorem 2.6 all hold with probability at least 1 — 0(n 



(5 > 2, if Co is replaced by C/3 for some numerical constant C > 0. 
3.2 Proof of Lemma l2"T8l 

We begin by introducing a piece of notation and set Zj = UV* — VrYj obeying 

Zj = (Ft - q- x VTP ai r T )Zj-x. 
Obviously Zj G T for all j > 0. First, note that when 

_ 2 pr log n 



> Cne" 



(for rectangular matrices, take q > Co e 2 jJrl ° g " (1) ), we have 

(2) 

II -^j" II 00 — ^W^j— 1|| 00 



(3.3) 



(3-4) 



by Lemma 3.1 (This holds with high probability because £lj and Zj—\ are independent, and this 
is why the golfing scheme is easy to use.) In particular, this gives that with high probability 



When q obeys the same estimate, 



\Zj\\oo < ei\\UV*\ 



by Theorem 2.6. In particular, this gives that with high probability 

\\Zj\\ F < e j \\UV*\\ F = e j Vr. 

Below, we will assume e < e" 1 . 



(3.5) 
(3.6) 
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Proof of (a). We prove the first part of the lemma and the argument parallels that in [22], see 
also [44]. From 

i 



we deduce 



\W L 



j 

j 



i-l 



<^(l-W^ll^lloo. 



The fourth step follows from Lemma 3.2 and the fifth from (3.5 ). Since ||C/V* || < ^/Jrr/n, this gives 

\\W L \\ < C'e 



for some numerical constant C whenever q obeys (3.3). 



Proof of (b). Since T n Y jo = 0, 

Vn(UV* + V T xY j0 ) = Vn(UV* - V T Y jo ) = Vn(Z jo ), 
and it follows from ( |3.6| ) that 

\\Z jo IIf < e jo \\UV*\\ F = e jo y/r. 
Since e < e _1 and jo > 21ogn, e j0 < 1/n 2 and this proves the claim. 

Proof of (c). We have UV* + W L = Zj + Yj and know that Yj is supported on J7 C . Therefore, 
since ||Zj ||f < A/8, it suffices to show that ||^' ||oo < A/8. We have 

||^ lloo<<r 1 ^||Pn j — 1 1| oo 

3 

< q 1 ll-^'-illoo 

3 

< g" 1 X) ^11^* Hoc- 
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Since ||C/V* ||oo < JJJr/n, this gives 



e 2 

iv- II < n' 

I -l JO II DO — ° " / ,, N 2 



yZ/ir (log n) 2 

for some numerical constant C whenever q obeys ( |3.3[ ). Since A = 1/y/n, [|3jf ||oo < A/8 if 

,/^r(logn) 2 \ 1/4 



e < C 



Summary. We have seen that (a) and (b) are satisfied if e is sufficiently small and jo > 2 log n. For 
(c), we can take e on the order of (pr(log n) 2 /n) 1 ^ , which will be sufficiently small as well provided 
that p r in (1.4) is sufficiently small. Note that everything is consistent since Cq e ~ 2 t irlo s n < ^ 



This concludes the proof of Lemma 2.8 



3.3 Proof of Lemma I2T91 

It is convenient to introduce the sign matrix E = sgn(So) distributed as 



1, w. p. p/2, 

0, w. p. 1 - p, (3.7) 
-1, w. p. p/2. 



We shall be interested in the event {H'Pq'Pt || < o"} which holds with large probability when a 



- v /p + e, see Corollary 2.7 In particular, for any a > 0, {||7'n7 ? T|| < c} holds with high probability 



provided p is sufficiently small. 

Proof of (a). By construction, 

W s = XV T ±E + XV T ± X)(7>nWn) fc £ 

k>l 

:=V T xW$ + V T xWf. 

For the first term, we have ||7 7 yxW , Q || < \\Wq \\ = A||-E||. Then standard arguments about the norm 
of a matrix with i.i.d. entries give [48] 

\\E\\ < ^^pnp 

with large probability. Since A = this gives ||W(f|| < 4y / p. When the matrix is rectangular, 

we have 

||£|| <4 V ^I^ 

with high probability. Since A = l/^/ri^ in this case, ||Wjf || < 4y / p as well. 

Set 1Z = J2k>i(.T~ > n'PT'Pn) k and observe that 1Z is self-adjoint. For the second term, \\V T ± Wf\\ < 
\\WfW, where Wf = X1Z(E). We need to bound the operator norm of the matrix 1Z(E), and use a 
standard covering argument to do this. Throughout, N denotes an 1/2-net for S n_1 of size at most 
6 n (such a net exists, see [30, Theorem 4.16]). Then a standard argument [48] shows that 

\\K(E)\\ = sup (y,TZ(E)x)<4 sup (y,TZ{E)x). 

1,2/eS" -1 x,y£N 
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For a fixed pair (x, y) of unit-normed vectors in TV x N, define the random variable 

X(x,y) := (y,K(E)x) = (K(yx*),E). 
Conditional on 0, = supp(-E'), the signs of E are i.i.d. symmetric and Hoeffding's inequality gives 

/ 2i 2 \ 

Now since ||?/x*||f = 1, the matrix lZ(yx*) obeys ||7^(yx*)||ir < \\1Z\\ and, therefore, 

P(sup \X(x,y)\>t\ri) <2|Af exp(-^A 

^x,yeN > V \\l<-\\ J 

Hence, 



P(||7e(i?)||>t|^)<2|iV| 2 exp(-^). 
On the event {\\VnV T \\ < a], 



L^t 1-CT 2 

fe>l 



and, therefore, unconditionally, 

P(||tt(£OII > «) < 2|iV| 2 exp(-^) + F(\\VnV T \\ >*), 7 = 

This gives 

F(\\\K(E)\\ >t)<2x 6 2 "exp(-^) +F(\\V n V T \\ > a). 

With A = 1/y/n, 

\\W S \\ < 1/4, 

with large probability, provided that a, or equivalently p, is small enough. 

Proof of (b). Observe that 

V n ±W s = -XV q ±Vt(Vq - VnP T Vnr l E. 
Now for (i, j) £ 0,°, = (ej, W s ej) = (eie*,W s ), and we have 

W§ = \(X(i,j),E), 

where X(i,j) is the matrix — (Vn — VnVT'Pn)~ 1 VnVT(eie*). Conditional on 17 = supp(E'), the 
signs of E are i.i.d. symmetric, and Hoeffding's inequality gives 

P(| W g|>tt|ll)<2«p(- ]p ^ jj5 ), 

and, thus, 

pfsup|Wg| > t\\n) < 2n 2 exp( 2 * A 
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Since (2.2) holds, we have 



\\PnPT{eie$\\ F < \\raRr\\\\r T {.eie^\\ F < a^+vrjn 
on the event {H'Pn'PT || < cr}. On the same event, \\(Vn - VnTrVn)" 1 \\ < (l-a 2 )- 1 and, therefore, 



2a 2 fir 



\X(iJ)\\ 2 F< 



{I -a 2 ) 2 n ' 



Then unconditionally, 



2 2 / 2 \ 2 

sup \Wg\ > tX) <2n 2 exp(- n ^ N ) + P(|| V Q V T \\ > v), 7 " (1 ~ ° ' 



fir J " la 1 



This proves the claim when fir < p' r n(logn) 1 and p' r is sufficiently small. 

4 Numerical Experiments and Applications 

In this section, we perform numerical experiments corroborating our main results and suggesting 
their many applications in image and video analysis. We first investigate Principal Component 
Pursuit's ability to correctly recover matrices of various rank from errors of various density. We then 
sketch applications in background modeling from video and removing shadows and specularities 
from face images. 

While the exact recovery guarantee provided by Theorem |1.1| is independent of the particular 
algorithm used to solve Principal Component Pursuit, its applicability to large scale problems 
depends on the availability of scalable algorithms for nonsmooth convex optimization. For the 
experiments in this section, we use the an augmented Lagrange multiplier algorithm introduced 
m [32, 51] In Section § we describe this algorithm in more detail, and explain why it is our 
algorithm of choice for sparse and low-rank separation. 

One important implementation detail in our approach is the choice of A. Our analysis identifies 
one choice, A = 1/ y / max(ni, ttq), which works well for incoherent matrices. In order to illustrate the 
theory, throughout this section we will always choose A = l/-^/max(ni, 712). For practical problems, 
however, it is often possible to improve performance by choosing A according to prior knowledge 
about the solution. For example, if we know that S is very sparse, increasing A will allow us to 
recover matrices L of larger rank. For practical problems, we recommend A = 1/ y / max(ni, ri2) as 
a good rule of thumb, which can then be adjusted slightly to obtain the best possible result. 

4.1 Exact recovery from varying fractions of error 



We first verify the correct recovery phenomenon of Theorem 1.1 on randomly generated problems. 
We consider square matrices of varying dimension n = 500, . . . , 3000. We generate a rank-r matrix 
Lq as a product Lq = XY* where X and Y are n x r matrices with entries independently sampled 
from a jV(0, 1/n) distribution. So is generated by choosing a support set f2 of size k uniformly at 
random, and setting So = VnE, where E is a matrix with independent Bernoulli ±1 entries. 

Table [l] (top) reports the results with r = rank(Lo) = 0.05 x n and k = ||5o||o = 0.05 x n 2 . 
Table [T] (bottom) reports the results for a more challenging scenario, rank(Lo) = 0.05 x n and 

8 Both [32, 51] have posted a version of their code online. 
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Dimension n 


rank(Lo) 


•So o 


rank(L) 


IpIIo 


II f T r, II n 

1 1 — -M) \\F 

ll-^o \\f 


# SVD 


Time(s) 


500 


25 


12,500 


25 


12,500 


1.1 x 10~ 6 


16 


2.9 


1,000 


50 


50,000 


50 


50,000 


1.2 x 10~ 6 


16 


12.4 


2,000 


100 


200,000 


100 


200,000 


1.2 x 10~ 6 


16 


61.8 


3,000 


250 


450,000 


250 


450,000 


2.3 x 10~ 6 


15 


185.2 


rank(Lo) = 0.05 x n, 


l-Sbllo = 0.05 x n 2 . 


Dimension n 


rank(Lo) 


l<sbl|o 


rank(L) 


IpIIo 


\\L-L \\ F 
ll-^o \\f 


# SVD 


Time(s) 


500 


25 


25,000 


25 


25,000 


1.2 x 10" 6 


17 


4.0 


1,000 


50 


100,000 


50 


100,000 


2.4 x 10~ 6 


16 


13.7 


2,000 


100 


400,000 


100 


400,000 


2.4 x 10~ 6 


16 


64.5 


3,000 


150 


900,000 


150 


900,000 


2.5 x 10~ 6 


16 


191.0 


rank(Lo) = 0.05 x n, 


|5o||o = 0.10 x n 2 . 



Table 1: Correct recovery for random problems of varying size. Here, Lq = XY* S R raxn 
with X,Y e R nxr ; X,Y have entries i.i.d. W(0, 1/n). S Q G {-1,0, 1}" X ™ has support chosen 
uniformly at random and independent random signs; ||£b||o is the number of nonzero entries 
in 5*o- Top: recovering matrices of rank 0.05 x n from 5% gross errors. Bottom: recovering 
matrices of rank 0.05 x n from 10% gross errors. In all cases, the rank of L and ^ - norrn of 
S'o are correctly estimated. Moreover, the number of partial singular value decompositions (# 
SVD) required to solve PCP is almost constant. 



k = 0.10 x n 2 . In all cases, we set A = 1 j \fn. Notice that in all cases, solving the convex PCP gives 
a result {L,S) with the correct rank and sparsity. Moreover, the relative error \\L — Lq\\f/\\Lq\\f 
is small, less than 10 -5 in all examples considered]^] 

The last two columns of Table [T] give the number of partial singular value decompositions 
computed in the course of the optimization SVD) as well as the total computation time. This 
experiment was performed in Matlab on a Mac Pro with dual quad-core 2.66 GHz Intel Xenon 
processors and 16 GB RAM. As we will discuss in Section [5] the dominant cost in solving the 
convex program comes from computing one partial SVD per iteration. Strikingly, in Table [T] the 
number of SVD computations is nearly constant regardless of dimension, and in all cases less than 
17, 10 This suggests that in addition to being theoretically well-founded, the recovery procedure 
advocated in this paper is also reasonably practical. 



4.2 Phase transition in rank and sparsity 

Theorem |1.1| shows that convex programming correctly recovers an incoherent low-rank matrix 
from a constant fraction p s of errors. We next empirically investigate the algorithm's ability to 
recover matrices of varying rank from errors of varying sparsity. We consider square matrices of 

9 We measure relative error in terms of L only, since in this paper we view the sparse and low-rank decomposition 
as recovering a low-rank matrix Lo from gross errors. So is of course also well-recovered: in this example, the relative 
error in S is actually smaller than that in L. 

10 One might reasonably ask whether this near constant number of iterations is due to the fact that random problems 
are in some sense well-conditioned. There is some validity to this concern, as we will see in our real data examples. [32] 
suggests a continuation strategy (there termed "Inexact ALM") that produces qualitatively similar solutions with a 
similarly small number of iterations. However, to the best of our knowledge its convergence is not guaranteed. 
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dimension n\ = ri2 = 400. We generate low-rank matrices Lq = XY* with X and Y independently 
chosen n x r matrices with i.i.d. Gaussian entries of mean zero and variance 1/n. For our first 
experiment, we assume a Bernoulli model for the support of the sparse term So, with random signs: 
each entry of So takes on value with probability 1 — p, and values ±1 each with probability p/2. 
For each (r, p) pair, we generate 10 random problems, each of which is solved via the algorithm of 
Sectionpl We declare a trial to be successful if the recovered L satisfies \\L — Lo\\f /\\Lo\\f < 10 -3 . 
Figure nT (left) plots the fraction of correct recoveries for each pair (r,p). Notice that there is a 
large region in which the recovery is exact. This highlights an interesting aspect of our result: the 
recovery is correct even though in some cases HSoIIf ^ 1 1 -^o I If (e.g., for r/n = p, \\Sq\\f is ^Jn = 20 



times larger!). This is to be expected from Lemma 2.4 the existence (or non-existence) of a dual 
certificate depends only on the signs and support of So and the orientation of the singular spaces 
of L . 

However, for incoherent Lq, our main result goes one step further and asserts that the signs of 
So are also not important: recovery can be guaranteed as long as its support is chosen uniformly 
at random. We verify this by again sampling Lq as a product of Gaussian matrices and choosing 
the support f2 according to the Bernoulli model, but this time setting So = 'PQSgn(Lo). One might 
expect such So to be more difficult to distinguish from Lq. Nevertheless, our analysis showed that 
the number of errors that can be corrected drops by at most 1/2 when moving to this more difficult 
model. Figure [T] (middle) plots the fraction of correct recoveries over 10 trials, again varying r and 
p. Interestingly, the region of correct recovery in Figure [T] (middle) actually appears to be broader 
than that in Figure [T] (left). Admittedly, the shape of the region in the upper-left corner is puzzling, 
but has been 'confirmed' by several distinct simulation experiments (using different solvers). 

Finally, inspired by the connection between matrix completion and robust PCA, we compare 
the breakdown point for the low-rank and sparse separation problem to the breakdown behavior of 
the nuclear- norm heuristic for matrix completion. By comparing the two heuristics, we can begin 
to answer the question how much is gained by knowing the location $7 of the corrupted entries'! 
Here, we again generate Lq as a product of Gaussian matrices. However, we now provide the 
algorithm with only an incomplete subset M = V^Lq of its entries. Each is included in O, 
independently with probability 1 — p, so rather than a probability of error, here, p stands for the 
probability that an entry is omitted. We solve the nuclear norm minimization problem 

minimize \\L\\^ subject to Vq±L = Vfi±M 

using an augmented Lagrange multiplier algorithm very similar to the one discussed in Section [5] 
We again declare Lq to be successfully recovered if \\L — Lo\\p/\\Lo\\f < 10~ 3 . Figure [l] (right) plots 
the fraction of correct recoveries for varying r, p. Notice that nuclear norm minimization successfully 
recovers Lq over a much wider range of (r,p). This is interesting because in the regime of large 
k, k = J7(n 2 ), the best performance guarantees for each heuristic agree in their order of growth 
- both guarantee correct recovery for rank(Lo) = 0(n/ log 2 n). Fully explaining the difference in 
performance between the two problems may require a sharper analysis of the breakdown behavior 
of each. 

4.3 Application sketch: background modeling from surveillance video 

Video is a natural candidate for low-rank modeling, due to the correlation between frames. One 
of the most basic algorithmic tasks in video surveillance is to estimate a good model for the 
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Figure 1: Correct recovery for varying rank and sparsity. Fraction of correct recoveries 
across 10 trials, as a function of rank(Lo) (x-axis) and sparsity of Sq (y-axis). Here, n\ = n% = 
400. In all cases, Lq — XY* is a product of independent n x r i.i.d. 7V(0, 1/n) matrices. Trials 
are considered successful if Lo||f/||£o||.f < 10~ 3 . Left: low-rank and sparse decomposition, 
sgn(5o) random. Middle: low-rank and sparse decomposition, Sq = VnSgii(L ) . Right: matrix 
completion. For matrix completion, p s is the probability that an entry is omitted from the 
observation. 



background variations in a scene. This task is complicated by the presence of foreground objects: 
in busy scenes, every frame may contain some anomaly. Moreover, the background model needs to 
be flexible enough to accommodate changes in the scene, for example due to varying illumination. 
In such situations, it is natural to model the background variations as approximately low rank. 
Foreground objects, such as cars or pedestrians, generally occupy only a fraction of the image 
pixels and hence can be treated as sparse errors. 

We investigate whether convex optimization can separate these sparse errors from the low- 
rank background. Here, it is important to note that the error support may not be well-modeled 
as Bernoulli: errors tend to be spatially coherent, and more complicated models such as Markov 
random fields may be more appropriate [11,52]. Hence, our theorems do not necessarily guarantee 
the algorithm will succeed with high probability. Nevertheless, as we will see, Principal Component 
Pursuit still gives visually appealing solutions to this practical low-rank and sparse separation 
problem, without using any additional information about the spatial structure of the error. 

We consider two example videos introduced in [31]. The first is a sequence of 200 grayscale 
frames taken in an airport. This video has a relatively static background, but significant foreground 
variations. The frames have resolution 176 x 144; we stack each frame as a column of our matrix 
M 6 u 25 > 344x200 . "We decompose M into a low-rank term and a sparse term by solving the convex 



PCP problem (1.1) with A = \j^fn{. On a desktop PC with a 2.33 GHz Core2 Duo processor 
and 2 GB RAM, our Matlab implementation requires 806 iterations, and roughly 43 minutes to 
converge Figure |2|a) shows three frames from the video; (b) and (c) show the corresponding 
columns of the low rank matrix L and sparse matrix S (its absolute value is shown here) . Notice 
that L correctly recovers the background, while S correctly identifies the moving pedestrians. The 
person appearing in the images in L does not move throughout the video. 



11 The paper [32] suggests a variant of ALM optimization procedure, there termed the "Inexact ALM" that finds 
a visually similar decomposition in far fewer iterations (less than 50). However, since the convergence guarantee for 
that variant is weak, we choose to present the slower, exact result here. 
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(a) Original frames (b) Low-rank L (c) Sparse S (d) Low-rank L (e) Sparse S 

Convex optimization (this work) Alternating minimization [47] 



Figure 2: Background modeling from video. Three frames from a 200 frame video sequence 
taken in an airport [31]. (a) Frames of original video M. (b)-(c) Low-rank L and sparse 
components S obtained by PCP, (d)-(e) competing approach based on alternating minimization 
of an rn-estimator [47]. PCP yields a much more appealing result despite using less prior 
knowledge. 



Figure[2](d) and (e) compares the result obtained by Principal Component Pursuit to a state-of- 
the-art technique from the computer vision literature, [47] . 12 That approach also aims at robustly 
recovering a good low-rank approximation, but uses a more complicated, nonconvex m-estimator, 
which incorporates a local scale estimate that implicitly exploits the spatial characteristics of natural 
images. This leads to a highly nonconvex optimization, which is solved locally via alternating 
minimization. Interestingly, despite using more prior information about the signal to be recovered, 
this approach does not perform as well as the convex programming heuristic: notice the large 
artifacts in the top and bottom rows of Figure [2] (d). 

In Figure [3j we consider 250 frames of a sequence with several drastic illumination changes. 
Here, the resolution is 168 x 120, and so M is a 20, 160 x 250 matrix. For simplicity, and to 



illustrate the theoretical results obtained above, we again choose A = X/^/ni, 13 For this example, 
on the same 2.66 GHz Core 2 Duo machine, the algorithm requires a total of 561 iterations and 36 
minutes to converge. 

Figure [3] (a) shows three frames taken from the original video, while (b) and (c) show the 
recovered low-rank and sparse components, respectively. Notice that the low-rank component 
correctly identifies the main illuminations as background, while the sparse part corresponds to the 



12 We use the code package downloaded from http://www.salleurl.edu/~ftorre/papers/rpca/rpca.zip, modi- 
fied to choose the rank of the approximation as suggested in [47] . 

13 For this example, slightly more appealing results can actually be obtained by choosing larger A (say, 2/a/ui). 
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(a) Original frames (b) Low-rank L (c) Sparse S (d) Low-rank L (e) Sparse S 

Convex optimization (this work) Alternating minimization [47] 

Figure 3: Background modeling from video. Three frames from a 250 frame sequence taken in 
a lobby, with varying illumination [31]. (a) Original video M. (b)-(c) Low-rank L and sparse S 
obtained by PCP. (d)-(e) Low-rank and sparse components obtained by a competing approach 
based on alternating minimization of an m-estimator [47]. Again, convex programming yields 
a more appealing result despite using less prior information. 

motion in the scene. On the other hand, the result produced by the algorithm of [47] treats some 
of the first illumination as foreground. PCP again outperforms the competing approach, despite 
using less prior information. These results suggest the potential power for convex programming as 
a tool for video analysis. 

Notice that the number of iterations for the real data is typically higher than that of the 
simulations with random matrices given in Table [TJ The reason for this discrepancy might be 
that the structures of real data could slightly deviate from the idealistic low-rank and sparse 
model. Nevertheless, it is important to realize that practical applications such as video surveillance 
often provide additional information about the signals of interest, e.g. the support of the sparse 
foreground is spatially piecewise contiguous, or even impose additional requirements, e.g. the 
recovered background needs to be non-negative etc. We note that the simplicity of our objective 
and solution suggests that one can easily incorporate additional constraints and more accurate 
models of the signals so as to obtain much more efficient and accurate solutions in the future. 

4.4 Application sketch: removing shadows and specularities from face images 

Face recognition is another problem domain in computer vision where low-dimensional linear models 
have received a great deal of attention. This is mostly due to the work of Basri and Jacobs, who 
showed that for convex, Lambertian objects, images taken under distant illumination lie near an 
approximately nine-dimensional linear subspace known as the harmonic plane [1]. However, since 
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(a) M 



(b) L 



(c) S 



(a) M 



(b) L 



(c) S 



Figure 4: Removing shadows, specularities, and saturations from face images, (a) Cropped 
and aligned images of a person's face under different illuminations from the Extended Yale 
B database. The size of each image is 192 x 168 pixels, a total of 58 different illuminations 
were used for each person, (b) Low-rank approximation L recovered by convex programming, 
(c) Sparse error S corresponding to specularities in the eyes, shadows around the nose region, 
or brightness saturations on the face. Notice in the bottom left that the sparse term also 
compensates for errors in image acquisition. 



faces are neither perfectly convex nor Lambertian, real face images often violate this low-rank 
model, due to cast shadows and specularities. These errors are large in magnitude, but sparse in 
the spatial domain. It is reasonable to believe that if we have enough images of the same face, 
Principal Component Pursuit will be able to remove these errors. As with the previous example, 
some caveats apply: the theoretical result suggests the performance should be good, but does not 
guarantee it, since again the error support does not follow a Bernoulli model. Nevertheless, as we 
will see, the results are visually striking. 

Figure [4] shows two examples with face images taken from the Yale B face database [18]. Here, 
each image has resolution 192 x 168; there are a total of 58 illuminations per subject, which we 
stack as the columns of our matrix M E K 32,256x58 . "We again solve PCP with A = 1/ \fn\. In 
this case, the algorithm requires 642 iterations to converge, and the total computation time on the 
same Core 2 Duo machine is 685 seconds. 

Figure [4] plots the low rank term L and the magnitude of the sparse term S obtained as the 
solution to the convex program. The sparse term S compensates for cast shadows and specular 
regions. In one example (bottom row of Figure [4] left), this term also compensates for errors in image 
acquisition. These results may be useful for conditioning the training data for face recognition, as 
well as face alignment and tracking under illumination variations. 
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5 Algorithms 



Theorem 1 1 . 1 1 shows that incoherent low- rank matrices can be recovered from nonvanishing fractions 
of gross errors in polynomial time. Moreover, as the experiments in the previous section attest, 
the low computation cost is guaranteed not only in theory, the efficiency is becoming practical for 
real imaging problems. This practicality is mainly due to the rapid recent progress in scalable 
algorithms for nonsmooth convex optimization, in particular for minimizing the l\ and nuclear 
norms. In this section, we briefly review this progress, and discuss our algorithm of choice for this 
problem. 

For small problem sizes, Principal Component Pursuit 

minimize + AHS^Ii 

subject to L + S = M 

can be performed using off-the-shelf tools such as interior point methods [21]. This was suggested 
for rank minimization in [16,45] and for low-rank and sparse decomposition [12] (see also [35]). 
However, despite their superior convergence rates, interior point methods are typically limited to 
small problems, say n < 100, due to the 0(n 6 ) complexity of computing a step direction. 

The limited scalability of interior point methods has inspired a recent flurry of work on first-order 
methods. Exploiting an analogy with iterative thresholding algorithms for ^i-minimization [49,50], 
Cai et. al. developed an algorithm that performs nuclear-norm minimization by repeatedly shrinking 
the singular values of an appropriate matrix, essentially reducing the complexity of each iteration 
to the cost of an SVD [6]. However, for our low-rank and sparse decomposition problem, this 
form of iterative thresholding converges slowly, requiring up to 10 4 iterations. Ma et. al. [20,36] 
suggest improving convergence using continuation techniques, and also demonstrate how Bregman 
iterations [41] can be applied to nuclear norm minimization. 

The convergence of iterative thresholding has also been greatly improved using ideas from 
Nesterov's optimal first-order algorithm for smooth minimization [37] , which was extended to non- 
smooth optimization in [2,38], and applied to ^i-minimization in [2,3,39]. Based on [2], Toh et. 
al. developed a proximal gradient algorithm for matrix completion which they termed Accelerated 
Proximal Gradient (APG). A very similar APG algorithm was suggested for low-rank and sparse 
decomposition in [33] . That algorithm inherits the optimal 0(1/ k 2 ) convergence rate for this class of 
problems. Empirical evidence suggests that these algorithms can solve the convex PCP problem at 
least 50 times faster than straightforward iterative thresholding (for more details and comparisons, 
see [33]). 

However, despite its good convergence guarantees, the practical performance of APG depends 
strongly on the design of good continuation schemes. Generic continuation does not guarantee good 



accuracy and convergence across a wide range of problem settings. In this paper, we have chosen 



to instead solve the convex PCP problem (1.1) using an augmented Lagrange multiplier (ALM) 
algorithm introduced in [32,51]. In our experience, ALM achieves much higher accuracy than 
APG, in fewer iterations. It works stably across a wide range of problem settings with no tuning 
of parameters. Moreover we observe an appealing (empirical) property: the rank of the iterates 
often remains bounded by rank(Lo) throughout the optimization, allowing them to be computed 
especially efficiently. APG, on the other hand, does not have this property. 



14 In our experience, the optimal choice may depend on the relative magnitudes of the L and S terms and the 
sparsity of the corruption. 
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The ALM method operates on the augmented Lagrangian 

l(L,S,Y) = \\L\\* + \\\S\\i + {Y,M-L-S) + ^\\M-L-S\\ 2 F . (5.1) 

A generic Lagrange multiplier algorithm [5] would solve PCP by repeatedly setting (L k ,Sk) = 
argmin^s l(L, S, Yfc), and then updating the Lagrange multiplier matrix via Yfe+i = + /u(M — 
Lk — Sk). 

For our low-rank and sparse decomposition problem, we can avoid having to solve a sequence of 
convex programs by recognizing that min^ l(L, S, Y) and mins l(L, S, Y) both have very simple and 
efficient solutions. Let S T : R — > R denote the shrinkage operator S T [x] = sgn(x) max(|x| — r, 0), 
and extend it to matrices by applying it to each element. It is easy to show that 

argmin/(L,S,y) = Sx^M — L + n~ l Y). (5.2) 

Similarly, for matrices X, let T> T {X) denote the singular value thresholding operator given by 
T> T (X) = US t {Tj)V* , where X = UT,V* is any singular value decomposition. It is not difficult to 
show that 

argmm/(L,S,y) = V lx (M-S-yT l Y). (5.3) 

Thus, a more practical strategy is to first minimize I with respect to L (fixing S), then minimize I 
with respect to S (fixing L), and then finally update the Lagrange multiplier matrix Y based on 
the residual M — L — S, a strategy that is summarized as Algorithm [T] below. 

Algorithm 1 (Principal Component Pursuit by Alternating Directions [32,51]) 

1: initialize: So = Yq = 0, /i > 0. 

2: while not converged do 
3: compute L k+1 = V^{M - S k - ^Yk); 
4: compute S k+ i = S\^(M - L k+ i + /U" 1 !^); 
5: compute Y k+1 =Y k + n(M - L k+1 - S k+1 ); 

6: end while 

7: output: L, S. 



Algorithm 1 is a special case of a more general class of augmented Lagrange multiplier algorithms 
known as alternating directions methods [51]. The convergence of these algorithms has been well- 
studied (see e.g. [29,34] and the many references therein, as well as discussion in [32,51]). Algorithm 
[T] performs excellently on a wide range of problems: as we saw in Section 3, relatively small 
numbers of iterations suffice to achieve good relative accuracy. The dominant cost of each iteration 
is computing Lk+i via singular value thresholding. This requires us to compute those singular 
vectors of M — Sk — [i~ x Yk whose corresponding singular values exceed the threshold \i. Empirically, 
we have observed that the number of such large singular values is often bounded by rank(Lo), 



allowing the next iterate to be computed efficiently via a partial SVD. 15 The most important 
implementation details for this algorithm are the choice of \x and the stopping criterion. In this 
work, we simply choose [i = nin2/4[|M||i, as suggested in [51]. We terminate the algorithm when 
\\M — L — S\\ F < S\\M\\ F , with 5 = 1(T 7 . 



15 Further performance gains might be possible by replacing this partial SVD with an approximate SVD, as suggested 
in [20] for nuclear norm minimization. 
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Very similar ideas can be used to develop simple and effective augmented Lagrange multiplier 



algorithms for matrix completion [32] , and for the robust matrix completion problem ( 1.5 1 discussed 



in Section 1.6 with similarly good performance. In the preceding section, all simulations and 
experiments are therefore conducted using ALM-based algorithms. For a more thorough discussion, 
implementation details and comparisons with other algorithms, please see [32,51]. 

6 Discussion 

This paper delivers some rather surprising news: one can disentangle the low-rank and sparse 
components exactly by convex programming, and this provably works under very broad conditions 
that are much broader than those provided by the best known results. Further, our analysis has 
revealed rather close relationships between matrix completion and matrix recovery (from sparse 
errors) and our results even generalize to the case when there are both incomplete and corrupted 



entries (i.e. Theorem 1.2). In addition, Principal Component Pursuit does not have any free 
parameter and can be solved by simple optimization algorithms with remarkable efficiency and 
accuracy. More importantly, our results may point to a very wide spectrum of new theoretical and 
algorithmic issues together with new practical applications that can now be studied systematically. 

Our study so far is limited to the low-rank component being exactly low-rank, and the sparse 
component being exactly sparse. It would be interesting to investigate when either or both these 
assumptions are relaxed. One way to think of this is via the new observation model M = Lq + So + 
iVo, where Nq is a dense, small perturbation accounting for the fact that the low-rank component 
is only approximately low-rank and that small errors can be added to all the entries (in some sense, 
this model unifies the classical PCA and the robust PCA by combining both sparse gross errors and 
dense small noise). The ideas developed in [7] in connection with the stability of matrix completion 
under small perturbations may be useful here. Even more generally, the problems of sparse signal 
recovery, low-rank matrix completion, classical PCA, and robust PCA can all be considered as 
special cases of a general measurement model of the form 

M = A(L )+B(S )+C(N ) 1 

where A, B, C are known linear maps. An ambitious goal might be to understand exactly under what 
conditions, one can effectively retrieve or decompose Lq and So from such noisy linear measurements 
via convex programming. 

The remarkable ability of convex optimizations in recovering low-rank matrices and sparse 
signals in high-dimensional spaces suggest that they will be a powerful tool for processing massive 
data sets that arise in image/video processing, web data analysis, and bioinformatics. Such data 
are often of millions or even billions of dimensions so the computational and memory cost can be 
far beyond that of a typical PC. Thus, one important direction for future investigation is to develop 
algorithms that have even better scalability, and can be easily implemented on the emerging parallel 
and distributed computing infrastructures. 
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7 Appendix 

7.1 Equivalence of sampling models 

We begin by arguing that a recovery result under the Bernoulli model automatically implies a 
corresponding result for the uniform model. Denote by Punif(m) an d PBer(p) probabilities calculated 
under the uniform and Bernoulli models and let "Success" be the event that the algorithm succeeds. 
We have 

n 2 

IPBer(p) (Success) = ^ P Ber(p ) (Success | |0| = k) P B er(p) (I ^1 = k) 
k=0 

rn—1 n 2 

< Yl P Bcr( P )(l^l = k) + >T IPunif(fc) (Success) P Bcr(p) (|0| = k) 
k=0 k=m 

< PBer(p)(|^| < m ) +IPunif(m) (Success), 

where we have used the fact that for k > m, PuniffJfe) (Success) < Punif(m) (Success), and that the 
conditional distribution of 0, given its cardinality is uniform. Thus, 

IPumf(m) (Success) > P B er( P ) (Success) -P Ber ( p )(|^| < m). 

z 2 n 2 

Take p = m/n + e, where e > 0. The conclusion follows from P Ber ( p )(|r2| < m) < e 2 p . In the 
other direction, the same reasoning gives 

rn 

IPBer(p) (Success) > ^ P Ber(p ) (Success | |fi| = fc)P Be r(p)(|fi| = k) 
k=0 

m 

> Punif(m) (Success) ^P Ber(p )( = k) 

k=0 

= Punif(m) (Success) P(|0 1 < 77l), 

and choosing m such that P(|^| > m) is exponentially small, establishes the claim. 

7.2 Proof of Lemma 13.11 

The proof is essentially an application of Bernstein's inequality, which states that for a sum of 
uniformly bounded independent random variables with |Y& — E Y&| < c, 

n t 2 
P (£(Y k - £ Y t )>t)<2ew(-^ T ^), (7.1) 

where a 2 is the sum of the variances, a 2 = Ylk=l Var(Yfc). 

Define via £Iq = ■ Sij = 1} where {<%} is an independent sequence of Bernoulli 

variables with parameter po- With this notation, Z' = Z — p„ Vt'PqqZ is given by 

Z' = E(1- Po%)^Me ie *) 
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so that Z'i j is a sum of independent random variables, 



Z 'iojo = 2 Yi i> Yi i = ( 1 ~ Po 1 5ij)Zij{VT(e i e* j ), e io e* jo ) 



We have 



£ Var(*S;) = (1 - po)po 1 £ l^fl^T^e*), e, e* )| 2 
< (I - Po )Po 1 \\Z\\lY,\^ e h V T(e to e* ))\ 2 



= (1- p )p~ 1 \\Zf oo \\V T (ei e* )\\ 2 F 

< (l-Po)Po \\ Z Woo — , 



< 



where the last inequality holds because of (2.2 ). Also, it follows from ( 1.2 ) that [(^(eie*), ej e* ) 
||^T( e i e j)IMI^r( e io e j )ll-F — 2pr/n so that \Yij\ < pQ 1 \\Z\\ OD pr /n. Then Bernstein's inequality 

P(|4|> e ||Z|| 0O )<2exp(-4^ 
J V lb pr 



J\\r II' i V-Wjo 

gives 



If po is as in Lemma 3.1 the union bound proves the claim. 
7.3 Proof of Theorem PI 

This section presents a proof of Theorem 1.2 which resembles that of Theorem |1.1[ Here and 
below, S' = Vn obs So so that the available data are of the form Y = Vn ohs Lo + S' . We make three 
observations. 

• If PCP correctly recovers Lq from the input data Vn obs Lo + S' (note that this means that 
L = Lq and S = S' ), then it must correctly recover Lq from Vq^Lq + Sq, where S' ' is 
a trimmed version of S' . The proof is identical to that of our elimination result, namely, 
Theorem 2.2 The derandomization argument then applies and it suffices to consider the case 
where the signs of S' Q are i.i.d. symmetric Bernoulli variables. 

• It is of course sufficient to prove the theorem when each entry in £l a b s is revealed with 
probability po := 0.1, i.e. when Q Q b s ~ Ber(po)- 

• We establish the theorem in the case where n\ = n<i = n as slight modifications would give 
the general case. 

Further, there are now three index sets of interest: 

• o b s are those locations where data are available. 

• r C S^obs are those locations where data are available and clean; that is, VrY = VtLq. 

• fi = Slobs \ r are those locations where data are available but totally unreliable. 

The matrix S'q is thus supported on Q. If £l Q h s ~ Ber(po); then by definition, £1 ~ Ber(po r )- 
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Dual certification. We begin with two lemmas concerning dual certification. 

Lemma 7.1 Assume \\Vt±Vt\\ < 1- Then (Lq,S' ) is the unique solution if there is a pair (W,F) 
obeying 

UV* + W = X(sgn(S' ) + F), 
with V T W = 0, \\W\\ < 1, V r ±F = and WF]]^ < 1. 

The proof is about the same as that of Lemma |2.4| and is discussed in very brief terms. The idea 
is to consider a feasible perturbation of the form (Lo + Hl, S' — Hs) obeying Vn ohs HL = Vn ohs Hg, 
and show that this increases the objective functional unless Hl = Hs = 0. Then a sequence of 



steps similar to that in the proof of Lemma 2.4 establishes 

\\L + Fill* + \\\S' - H s \\i > \\L \U + A||^||i + (1 - P)(\\V T ±H L \\* + A||^F L ||i) 



(7.2) 



where /3 = max(||W||, llFlloo). Finally, \\P t ±Hl\\*+^\\PtHl\\i vanishes if and only Hl ET ± PiT = 
{0}. 

Lemma 7.2 Assume that for any matrix M, \\VTT' r ±M\\ F < ^ 1 1 "P^j- 'Pp^ iW 1 1 ^ and take X > 4/n. 
Then (Lq,S' ) is the unique solution if there is a pair (W,F) obeying 

UV* + W + V T D = X(sgn(S' ) + F), 

with V T W = 0, ||W|| < 1/2, V T ±F = and {{F^ < 1/2, and \\V T D\\ F < n~ 2 . 

Note that \\V T V T ±M\\ F < n\\V T ±V r ±M\\ F implies T 1 - n T = {0}, or equivalently \\V T ±V T \\ < 1- 
Indeed if M E T 1 - n T, V T V r ±M = M while V T ±V r ±M = 0, and thus M = 0. 



Proof It follows from (7.2) together with the same argument as in the proof of Lemma 7.2 that 



\Lq + Fi||* + X\\S' - H s \\i > \\Lq\U + A||^||i + J (\\V T ±H L \\* + X\\VrH L \\i) - ^\\V T n L \\ 



n- 



Observe now that 



\\VtH l \\f < WVtVtHlWf + \\PtV t ±H l \\f 

< WVtVtHlWf + n\\V T ^V T ^H L \\ F 

< WVtVtHlWf + n{\\V T ^TH L \\F + \\V T xH L \\ F ) 

< (n + l)\\VrH L \\ F + n\\P T ±H L \\ F . 



Using both \\VrH L \\ F < \\VtHl\\i and \\V t ±H l \\f < \\V T ±H L \U, we obtain 



\\L + F L ||* + X\\S' -H s \\i> \\L \U + A||^||i + 
The claim follows from T 1 - n T = {0}. 



1 1 



n 



\V t ±Hl\\* + 



A n + 1 



n- 



YPvHlWi. 



Lemma 7.3 Under the assumptions of Theorem 1.2, the assumption of Lemma 7.2 is satisfied with 



high probability. That is, ||'Pi'Pr±M||^ < n\\V F ±'P^±M\\F for all M. 
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Proof Set po = p (l - r) and M' = V r ±M. Since T ~ Ber(p ), Theorem |2~6l gives \\V T - 
Pq X VtPtPt\ < 1/2 with high probability. Further, because \\VtV t M'\\ f = \\VvV T ±M'\\ F , we 
have 

\\VtV t M'\\ f < \\P T ±M'\\ F . 

In the other direction, 

Pq X \\PtV t M'\\ 2 f = Pv 1 (VtM',VtVtV t M') 

= (V T M',V T M') + (V T M', (p^VtVtVt - V T )M') 

> \\V T M'\\ 2 F - h\\V T M'f F = \\\P T M'\\ 2 F . 

In conclusion, \\T> T ±M'\\ F > \\VvVtM' \\p > ®-\\V T M'\\ F , and the claim follows since f > ±. ■ 

Thus far, our analysis shows that to establish our theorem, it suffices to construct a pair 
(Y L , W s ) obeying 



\\V T ±Y L \\ < 1/4, 
\\V T Y L - UV*\\ F < n~ 2 , 
V r ±Y L = 0, 
||W L ||oo < A/4, 



and 



V T W S = 0, 
\\W S \\ < 1/4, 
V n W s = Asgn(^), 
V Q x W s = 0, 

obs 

ll|7W 5 ||oc < A/4. 



(7.3) 



Indeed, by definition, Y L + W s obeys 

Y L + W S = A(sgn(5^) + F), 



where F is as in Lemma 7.2 and it can also be expressed as 

Y L + W s = UV* + W + V T D, 
where W and VtD are as in this lemma as well. 

Construction of the dual certificate Y L . We use the golfing scheme to construct Y L . Think 
of T ~ Ber(po) with po = po(l — r) as ^i<j<j ^j, where the sets Tj ~ Ber(g) are independent, and 
q obeys po = 1 — (1 — q) J0 . Here, we take jo = [~31ogn] , and observe that q > po/jo as before. Then 
starting with Y$ = 0, inductively define 



and set 



Y 3 = 3^-1 + q-^VriUV* - Y^), 
Y L = Y j0 = ^ J>r^,-i, z i = ( p T ~ q^VrVr^Zj-L 



(7.4) 



By construction, V r ±Y L = 0. Now just as in Section (3.2), because q is sufficiently large, \\Zj\\ < 
e~ :, \\UV*\\ 00 and < e~ 3 yfr, both inequality holding with large probability. The proof is now 

identical to that in (2.5). First, the same steps show that 



V T ±Y L \\ < C\ 



' n log n , 



UV*\ 



c 



, //ir(logn) 2 
np 
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Whenever po > Co ^^ 1 ^"^ for a sufficiently large value of the constant Cq (which is possible 
provided that p r in (1.6) is sufficiently small), this terms obeys HT't-'-^'II < 1/4 as required. 
Second, 

\\V T Y L - UV*\\ F = \\Z j0 \\ F < e- 31ogn ^ < n- 2 . 
And third, the same steps give 



|^ L ||oo <g _1 ||^*||ooX) e-j <3(l-e" 1 ). 



' //r(logn) 5 
pin 2 



Now it suffices to bound the right-hand side by j = This is automatic when pq > 

/w(log n) 2 



Co „ whenever Cq is sufficiently large and, thus, the situation is as before. In conclusion, we 
have established that Y obeys (7.3) with high probability. 

Construction of the dual certificate W s . We first establish that with high probability, 

WVrVnW < V?Po, t' = t + t , (7.5) 

where To(r) is a continuous function of r approaching zero when r approaches zero. In other words, 
the parameter r' may become arbitrary small constant by selecting r small enough. This claim is 
a straight application of Corollary |2.7| We also have 



WP(iP (T+a xjV n \\ < 2t'. (7.6) 
with high probability. This second claim uses the identity 

V 1 obs-' 

This is well defined since the restriction of VtVu^Vt to T is invertible. Indeed, Theorem 2.6 gives 
VtV^Vt > ^V T and, therefore, \\{V T V aohs VT)- l \\ < 2p^. Hence, 

\\VnV {T+n ^ s) Vn\\ < 2 Po 1 \\Vnr T \\ 2 , 



and (7.6) follows from (7.5). 

Setting E = sgn(5 ), this allows to define W s via 

W s = A(J _ v (T+nLs) )(Vn - VnV {T+ ^ bs) 'Pn)- 1 E 

■■= (i-v^^jiwi + wf), 

where W$ = XE, and Wf = TIE with K = ^ k>1 (PnF \ T+n ± )V n ) k - The operator K is self- 
adjoint and obeys < 1 2t 2 —, with high probability. By construction, VtW s = V n ^ W s = 

obs 

and VqW s = Asgn(S' ). It remains to check that both events \\W S \\ < 1/4 and HPrW^loo < A/4 
hold with high probability. 

Control of \\W S \\. For the first term, we have \\(2 — V^ T+n ± ))W(f || < \\Wq\\ = X\\E\\. Because 
the entries of E are i.i.d. and take the value ±1 each with probability pqt/2, and the value with 
probability 1 — pqt, standard arguments give 

\\E\\ < 4 v / np (r + r ) 



35 



with large probability. Since A = I/^/pqu, \\Wq\\ < 4-^r + tq < 1/8 with high probability, provided 
r is small enough. 

For the second term, 'P( T+n ± -j)Wf\\ < X\\TZE\\, and the same covering argument as before 
gives 

F(\\\K(E)\\ >t)<2x6 2 " exp (-^) + P(||ft|| > a). 

Since A = 1 / yjnpo this shows that < 1/4 with high probability, since one can always choose 

a, or equivalently t' = r + To, sufficiently small. 
Control of ||7W S ||oo. For G T, we have 

^ = (e ie *,^ 5 ) = A(X(i,i),^), 

where 

X(i,j) = {Vn-VnV {T+nLs) Vnr l VnV {T+uLsV eie*. 
The same strategy as before gives 

P( sup |Wg| > <2n 2 exp(--M +p( sup j) || F > a) . 

It remains to control the Frobenius norm of X(i,j). To do this, we use the identity 

which gives 

WVnV^n^e^y < ^\\V T e t e*\\ F < 

with high probability. This follows from the fact that WiVrVn^Vr^W < 2y»Q 1 and WPnV T \\ < 
y/pQT 1 as we have already seen. Since we also have \\(Vn — 1- > n'P(T+n\ )'Pn)~ 1 \\ < \\ T ' m § n 
probability, 

sup ||X(z,j)|| F <-i— 
This shows that H^rW^Hoo < A/4 if t', or equivalently r, is sufficiently small. 
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